Mark Meyer | Parametric surfaces | colored trefoil knot
size(550, 550)
background(0)
nofill()
stroke(0.5, 0.25)
strokewidth(0.5)
translate(WIDTH/2, HEIGHT/2)
# n defines how many vertices we find: the number will be n^2.
# Large values of n can take a while to render!
n = 150
# --- TREFOIL EQUATION ------------------------------------------------------
# Define parametric equations for each axis.
# Although we define an z coordinate,
# it is only used in transformations and shading;
# it is not used in the drawing
# (although it could be if you wanted to add perspective).
r = 120
def x(u, v):
return r * sin(3*u) / (2 + cos(v))
def y(u, v):
return r * (sin(u) + 2 * sin(2*u)) / (2 + cos(v + 2*pi/3))
def z(u, v):
return r/2 * (cos(u) - 2 * cos(2*u)) \
* (2 + cos(v)) \
* (2 + cos(v + 2*pi/3)) / 4
def fit_to_domain(u, v):
u = -pi + u * 2*pi / n
v = -pi + v * 2*pi / (n-1)
return u, v
from Numeric import *
def matrix(u, v, index):
u, v = fit_to_domain(u, v)
return where(index==0, x(u, v),
where(index==1, y(u,v),
where(index==2, z(u,v), 1)
)
)
# --- MATRIX ROTATION ------------------------------------------------------
def rotate_matrix(matrix, x, y, z):
rot_x = array(
([1, 0, 0],
[0, cos(x), -sin(x)],
[0, sin(x), cos(x)]
), Float
)
rot_y = array(
([cos(y), 0, sin(y)],
[0, 1, 0],
[-sin(y), 0, cos(y)]
), Float
)
rot_z = array(
([cos(z), -sin(z), 0],
[sin(z), cos(z), 0],
[0, 0, 1]
), Float
)
matrix = matrixmultiply(matrix, rot_x)
matrix = matrixmultiply(matrix, rot_y)
matrix = matrixmultiply(matrix, rot_z)
return matrix
# --- VECTOR NORMALIZATION -------------------------------------------------
def length(vector):
return sqrt(vector[0]**2 + vector[1]**2 + vector[2]**2)
def normalize_vector(face):
""" Calculate the normal vector between -1 and 1 for a flat surface """
norm = cross_product((face[0,0] - face[1,1]), (face[1,0] - face[0,1]))
return norm / length(norm)
# --- ORTHOGRAPHIC PROJECTION ----------------------------------------------
def project(rows):
"""
Go through the array and draw rectangles.
There is probably a more efficent way to do this.
"""
for i in range(len(rows) -1):
for j in range(len(rows[i])-1):
face = rows[i:i+2, j:j+2]
normal = normalize_vector(face)
light_angle = (dot(normal, light))
if (light_angle < 0):
light_angle = 0
# Coloring scheme based on t between -1.0 and 1.0.
colormode(HSB)
t = float(i) / len(rows) * 2 - 1
fill(
0.6 + 0.2 * abs(t**2),
0.8 - 0.4 * light_angle,
0.3 + 0.7 * light_angle,
0.75
)
# Draw only forward facing facets.
if (dot(view, normal) > 0):
beginpath(
face[0, 0, 0],
face[0, 0, 1]
)
lineto(
face[1, 0, 0],
face[1, 0, 1]
)
lineto(
face[1, 1, 0],
face[1, 1, 1]
)
lineto(
face[0, 1, 0],
face[0, 1, 1]
)
path = endpath()
# Save z with path for later z sorting.
path.z = face[1,0,2]
# --- ROTATION SLIDERS -----------------------------------------------------
var("rot_x", NUMBER, 0.00, -pi, pi)
var("rot_y", NUMBER, pi/2, -pi, pi)
var("rot_z", NUMBER, 0.00, -pi, pi)
# --- LIGHTING SLIDERS -----------------------------------------------------
var("light_x", NUMBER, 1.0, -1.0, 1.0)
var("light_y", NUMBER, 1.0, -1.0, 1.0)
var("light_z", NUMBER, 1.0, -1.0, 1.0)
light = array([light_x, light_y, light_z], Float)
# --- FACET CULLING --------------------------------------------------------
# Draw only forward facing facets.
view = array([0, 0, 1], Float)
sphere = fromfunction(matrix, (n+1, n+1, 3))
sphere = rotate_matrix(sphere, rot_x, rot_y, rot_z)
project(sphere)
# z-sort to avoid problems with overlapping surfaces:
sorted_grobs = list(canvas)
sorted_grobs.sort(lambda p1, p2: cmp(p1.z, p2.z))
canvas._grobs = sorted_grobs
